Overview

Setup

(hide)

Click on tabs to display additional information.

Libraries

# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS)  # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(Maaslin2)
library(tidyverse)
library(RColorBrewer)  # For color palettes in heatmaps
library(pheatmap)      # For heatmap visualization

Plotting

# Define treatment order and color palette
treatment_order <- c(
  "A- T- P-",  # Control
  "A- T- P+",  # Parasite
  "A+ T- P-",  # Antibiotics
  "A+ T- P+",  # Antibiotics_Parasite
  "A- T+ P-",  # Temperature
  "A- T+ P+",  # Temperature_Parasite
  "A+ T+ P-",  # Antibiotics_Temperature
  "A+ T+ P+"   # Antibiotics_Temperature_Parasite
)

# Custom color palette matching treatment order
treatment_colors <- c(
  "#1B9E77",  # A- T- P- (Control)
  "#D95F02",  # A- T- P+ (Parasite)
  "#7570B3",  # A+ T- P- (Antibiotics)
  "#E7298A",  # A+ T- P+ (Antibiotics_Parasite)
  "#66A61E",  # A- T+ P- (Temperature)
  "#E6AB02",  # A- T+ P+ (Temperature_Parasite)
  "#A6761D",  # A+ T+ P- (Antibiotics_Temperature)
  "#666666"   # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)

# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)

Functions

# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
  samdat <- phyloseq::sample_data(ps)
  df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
  return(df)
}

# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
  ps <- microViz::ps_get(ps)
  df <- samdatAsDataframe(ps)
  df <- dplyr::rename(.data = df, ...)
  phyloseq::sample_data(ps) <- df
  return(ps)
}

# SourceFolder function
source(here::here("Code", "R", "Functions", "StartFunctions", "sourceFolder.R"))

# Import all helper functions found in `/Functions`
sourceFolder(here::here("Code", "R", "Functions", "StartFunctions"), T)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
## The following object is masked from 'package:coin':
## 
##     pvalue
## Loading required package: future
## 
## Attaching package: 'future'
## The following object is masked from 'package:survival':
## 
##     cluster
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:microViz':
## 
##     stat_chull
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
##    in the original pheatmap() are identically supported in the new function. You 
##    can still use the original function by explicitly calling pheatmap::pheatmap().
## 
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
## 
##     pheatmap
## Loading required package: data.table
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following object is masked from 'package:purrr':
## 
##     transpose
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## Loading required package: doParallel
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: iterators
## Loading required package: GUniFrac
## Registered S3 method overwritten by 'rmutil':
##   method         from
##   print.response httr
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: rlang
## 
## Attaching package: 'rlang'
## The following object is masked from 'package:magrittr':
## 
##     set_names
## The following object is masked from 'package:data.table':
## 
##     :=
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## Attaching package: 'microbiome'
## The following object is masked from 'package:vegan':
## 
##     diversity
## The following object is masked from 'package:scales':
## 
##     alpha
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:base':
## 
##     transform
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following object is masked from 'package:ggpubr':
## 
##     rotate
## The following object is masked from 'package:dplyr':
## 
##     where
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## 3 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/StartFunctions"
sourceFolder(here::here("Code", "R", "Functions", "HelperFunctions"), T)
## 9 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/HelperFunctions"
# sourceFolder(here::here("Code", "R", "Functions", "AnalysisScripts"), T)


# Maaslin2

# Define a modified version of run_maaslin2
run_maaslin2_modified <- function(tmp.PS, 
                         tmp.sample.names = "Sample",
                         tmp.rank = "Genus",
                         tmp.fixed,
                         tmp.reference,
                         tmp.output,
                         tmp.plot.scatter = T,
                         tmp.plot.heatmap = T,
                         tr = 'LOG', 
                         norm = 'TSS'){
  
  tmp.input.otu = tmp.PS %>% microViz::tax_agg(rank = "Genus") %>% otu_get() %>% as.data.frame() %>%
      setNames(gsub(" ", "_", names(.)))
  
  tmp.input.meta = microViz::samdat_tbl(tmp.PS) %>% as.data.frame()
  row.names(tmp.input.meta) <- tmp.input.meta[[tmp.sample.names]]
  
  # Create output directory if it doesn't exist
  dir.create(tmp.output, recursive = TRUE, showWarnings = FALSE)
  
  obj <- Maaslin2::Maaslin2(input_data = tmp.input.otu,
                            input_metadata = tmp.input.meta,
                            output = tmp.output, 
                            min_prevalence = 0,
                            normalization = norm, 
                            transform = tr,
                            fixed_effects = tmp.fixed, 
                            reference = tmp.reference,
                            standardize = F,
                            plot_heatmap = tmp.plot.heatmap, 
                            plot_scatter = tmp.plot.scatter,
                            cores = 8)
  
  # Process results with corrected column selection
  res <- obj$results %>%
    filter(metadata == 'group' & value == 'case') %>%
    mutate(df = N - length(tmp.fixed) - 1,
           method = paste0('MaAsLin2 (', tr, ', ', norm, ')')) %>%
    dplyr::select(feature, coef, stderr, df, pval, method) %>%
    rename(taxon = feature,
           est = coef,
           se = stderr,
           p = pval)
  
  return(res)
}

#### GT TABLE STYLES

apply_GT_tabStyles <- function(gt_table, columns, coef_columns) {
  for (i in seq_along(columns)) {
    column <- columns[i]
    coef_column <- coef_columns[i]
    
    coef_vals <- gt_table$`_data`[[coef_column]]
    colors <- gradient_palette(coef_vals)
    
    for (j in seq_along(colors)) {
      cell_color <- ifelse(is.na(coef_vals[j]), "grey95", colors[j])
      cell_text <- "black" # ifelse(is.na(coef_vals[j]), "black", "white")
      cell_size <- "medium"
      bold_text <- ifelse(gt_table$`_data`[[column]][j] %in% c("+", "-"), "bold", "normal")
      
      gt_table <- gt_table %>%
        tab_style(
          style = list(cell_fill(color = cell_color),
                       cell_text(color = cell_text,
                                 weight = bold_text,
                                 align = "center",
                                 size = cell_size
                       )
          ),
          locations = cells_body(
            columns = !!sym(column),
            rows = j
          )
        )
    }
  }
  gt_table
}

# Function to create treatment-specific heatmap from MaAsLin2 results
create_treatment_heatmap <- function(maaslin_results, top_n = 50, title = "Treatment-specific Taxa Abundance Changes") {
  
  # Set seed for reproducibility
  set.seed(42)
  
  # Filter for significant results and get top taxa
  heatmap_data <- maaslin_results %>%
    dplyr::filter(qval < 0.05) %>%
    dplyr::arrange(qval) %>%
    dplyr::slice_head(n = top_n) %>%
    dplyr::select(Taxon, value, coef, Phylum) %>%
    dplyr::distinct(Taxon, value, .keep_all = TRUE)  # Remove duplicates if any
  
  # Pivot to wide format with treatments as columns
  heatmap_matrix <- heatmap_data %>%
    tidyr::pivot_wider(
      id_cols = c(Taxon, Phylum),
      names_from = value,
      values_from = coef,
      values_fill = 0
    ) %>%
    dplyr::arrange(Phylum, Taxon) %>%
    tibble::column_to_rownames("Taxon")
  
  # Extract phylum information for row annotations
  phylum_info <- heatmap_matrix$Phylum
  heatmap_matrix <- heatmap_matrix %>%
    dplyr::select(-Phylum) %>%
    as.matrix()
  
  # Order columns according to treatment_order (only include treatments that exist in the data)
  existing_treatments <- colnames(heatmap_matrix)
  ordered_treatments <- treatment_order[treatment_order %in% existing_treatments]
  heatmap_matrix <- heatmap_matrix[, ordered_treatments, drop = FALSE]
  
  # Create row annotation data frame
  annotation_row <- data.frame(
    Phylum = phylum_info,
    row.names = rownames(heatmap_matrix)
  )
  
  # Create color palette for phyla
  phylum_colors <- setNames(
    RColorBrewer::brewer.pal(length(unique(phylum_info)), "Set3"),
    unique(phylum_info)
  )
  
  # Define annotation colors
  ann_colors <- list(
    Phylum = phylum_colors
  )
  
  # Create color palette for heatmap
  # Use a diverging color palette centered at 0
  max_abs_coef <- max(abs(heatmap_matrix), na.rm = TRUE)
  heatmap_colors <- colorRampPalette(c("blue", "white", "red"))(100)
  
  # Create the heatmap
  pheatmap::pheatmap(
    heatmap_matrix,
    color = heatmap_colors,
    breaks = seq(-max_abs_coef, max_abs_coef, length.out = 101),
    annotation_row = annotation_row,
    annotation_colors = ann_colors,
    cluster_rows = TRUE,
    cluster_cols = FALSE,  # Keep treatment order as is
    show_rownames = TRUE,
    show_colnames = TRUE,
    fontsize_row = 8,
    fontsize_col = 10,
    angle_col = 45,
    main = title,
    treeheight_row = 20,
    treeheight_col = 0
  )
}

# Function to create comparison tables for differentially abundant taxa
create_taxa_comparison_tables <- function(maaslin_results, section_title = "Taxa Comparison") {
  
  # Set seed for reproducibility
  set.seed(42)
  
  if (is.null(maaslin_results) || nrow(maaslin_results) == 0) {
    return(list(
      summary_table = NULL,
      overlap_table = NULL,
      direction_table = NULL,
      top10_table = NULL
    ))
  }
  
  # Filter for significant results only
  sig_results <- maaslin_results %>%
    dplyr::filter(qval < 0.05) %>%
    dplyr::select(Taxon, value, coef, Phylum, qval)
  
  if (nrow(sig_results) == 0) {
    return(list(
      summary_table = NULL,
      overlap_table = NULL,
      direction_table = NULL,
      top10_table = NULL
    ))
  }
  
  # 1. Summary table - Number of significant taxa per treatment
  summary_table <- sig_results %>%
    dplyr::group_by(value) %>%
    dplyr::summarise(
      Total_Significant = dplyr::n(),
      Positive_Effect = sum(coef > 0),
      Negative_Effect = sum(coef < 0),
      .groups = "drop"
    ) %>%
    dplyr::arrange(factor(value, levels = treatment_order)) %>%
    gt::gt() %>%
    gt::tab_header(
      title = paste("Summary of Significant Taxa by Treatment"),
      subtitle = section_title
    ) %>%
    gt::cols_label(
      value = "Treatment",
      Total_Significant = "Total Significant",
      Positive_Effect = "Positive Effect",
      Negative_Effect = "Negative Effect"
    ) %>%
    gt::tab_style(
      style = cell_text(weight = "bold"),
      locations = cells_column_labels()
    ) %>%
    gt::fmt_number(
      columns = c(Total_Significant, Positive_Effect, Negative_Effect),
      decimals = 0
    )
  
  # Color the Treatment column rows by treatment color
  for (i in seq_along(treatment_order)) {
    trt <- treatment_order[i]
    if (trt %in% unique(sig_results$value)) {
      trt_color <- treatment_color_scale[trt]
      summary_table <- summary_table %>%
        gt::tab_style(
          style = cell_fill(color = trt_color),
          locations = cells_body(
            columns = "value",
            rows = value == trt
          )
        ) %>%
        gt::tab_style(
          style = cell_text(color = "white", weight = "bold"),
          locations = cells_body(
            columns = "value",
            rows = value == trt
          )
        )
    }
  }
  
  # 2. Top 10 taxa table - Most significant taxa by treatment
  top10_data <- sig_results %>%
    dplyr::group_by(value) %>%
    dplyr::arrange(qval, .by_group = TRUE) %>%
    dplyr::slice_head(n = 10) %>%
    dplyr::ungroup() %>%
    dplyr::select(Taxon, value, coef, Phylum, qval) %>%
    dplyr::arrange(factor(value, levels = treatment_order), qval)
  
  if (nrow(top10_data) > 0) {
    top10_table <- top10_data %>%
      gt::gt() %>%
      gt::tab_header(
        title = "Top 10 Most Significant Taxa by Treatment",
        subtitle = paste(section_title, "- Ranked by q-value")
      ) %>%
      gt::cols_label(
        value = "Treatment",
        coef = "Coefficient",
        qval = "q-value"
      ) %>%
      gt::tab_style(
        style = cell_text(weight = "bold"),
        locations = cells_column_labels()
      ) %>%
      gt::fmt_scientific(
        columns = qval,
        decimals = 3
      ) %>%
      gt::fmt_number(
        columns = coef,
        decimals = 3
      ) %>%
      gt::tab_style(
        style = cell_text(weight = "bold"),
        locations = cells_body(columns = "Taxon")
      )
    
    # Color treatment rows by treatment color
    for (i in seq_along(treatment_order)) {
      trt <- treatment_order[i]
      if (trt %in% unique(top10_data$value)) {
        trt_color <- treatment_color_scale[trt]
        top10_table <- top10_table %>%
          gt::tab_style(
            style = cell_fill(color = trt_color),
            locations = cells_body(
              columns = "value",
              rows = value == trt
            )
          ) %>%
          gt::tab_style(
            style = cell_text(color = "white", weight = "bold"),
            locations = cells_body(
              columns = "value",
              rows = value == trt
            )
          )
      }
    }
  } else {
    top10_table <- NULL
  }
  
  # 3. Overlap table - Compare taxa between treatments
  treatments <- unique(sig_results$value)
  if (length(treatments) > 1) {
    # Create list of taxa for each treatment
    treatment_taxa <- list()
    for (trt in treatments) {
      treatment_taxa[[trt]] <- sig_results %>%
        dplyr::filter(value == trt) %>%
        dplyr::pull(Taxon)
    }
    
    # Create overlap matrix
    overlap_matrix <- matrix(0, nrow = length(treatments), ncol = length(treatments))
    rownames(overlap_matrix) <- treatments
    colnames(overlap_matrix) <- treatments
    
    for (i in seq_along(treatments)) {
      for (j in seq_along(treatments)) {
        if (i == j) {
          overlap_matrix[i, j] <- length(treatment_taxa[[treatments[i]]])
        } else {
          overlap_matrix[i, j] <- length(intersect(treatment_taxa[[treatments[i]]], 
                                                  treatment_taxa[[treatments[j]]]))
        }
      }
    }
    
    # Convert to data frame for GT table
    overlap_df <- as.data.frame(overlap_matrix) %>%
      tibble::rownames_to_column("Treatment") %>%
      dplyr::arrange(factor(Treatment, levels = treatment_order))
    
    overlap_table <- overlap_df %>%
      gt::gt() %>%
      gt::tab_header(
        title = "Taxa Overlap Between Treatments",
        subtitle = paste(section_title, "- Number of shared significant taxa")
      ) %>%
      gt::tab_style(
        style = cell_text(weight = "bold"),
        locations = cells_column_labels()
      ) %>%
      gt::fmt_number(
        columns = -Treatment,
        decimals = 0
      ) 
    
    # Apply treatment colors to column headers and row headers
    for (i in seq_along(treatments)) {
      trt <- treatments[i]
      if (trt %in% names(treatment_color_scale)) {
        trt_color <- treatment_color_scale[trt]
        
        # Color column header
        overlap_table <- overlap_table %>%
          gt::tab_style(
            style = cell_fill(color = trt_color),
            locations = cells_column_labels(columns = trt)
          ) %>%
          gt::tab_style(
            style = cell_text(color = "white", weight = "bold"),
            locations = cells_column_labels(columns = trt)
          )
        
        # Color row header (first column) for matching treatment
        overlap_table <- overlap_table %>%
          gt::tab_style(
            style = cell_fill(color = trt_color),
            locations = cells_body(
              columns = "Treatment",
              rows = Treatment == trt
            )
          ) %>%
          gt::tab_style(
            style = cell_text(color = "white", weight = "bold"),
            locations = cells_body(
              columns = "Treatment",
              rows = Treatment == trt
            )
          )
      }
    }
  } else {
    overlap_table <- NULL
  }
  
  # 4. Direction comparison table - Compare positive/negative effects
  if (length(treatments) > 1) {
    direction_data <- sig_results %>%
      dplyr::group_by(value) %>%
      dplyr::summarise(
        Positive_Taxa = list(Taxon[coef > 0]),
        Negative_Taxa = list(Taxon[coef < 0]),
        .groups = "drop"
      )
    
    # Create direction comparison matrix
    direction_matrix <- matrix("", nrow = length(treatments), ncol = length(treatments))
    rownames(direction_matrix) <- treatments
    colnames(direction_matrix) <- treatments
    
    for (i in seq_along(treatments)) {
      for (j in seq_along(treatments)) {
        if (i == j) {
          # For diagonal cells, show total positive and negative counts for that treatment
          # Find the correct treatment data
          trt_name <- treatments[i]
          trt_data <- direction_data[direction_data$value == trt_name, ]
          
          if (nrow(trt_data) > 0) {
            pos_count <- length(trt_data$Positive_Taxa[[1]])
            neg_count <- length(trt_data$Negative_Taxa[[1]])
            direction_matrix[i, j] <- paste0(
              "Pos: ", pos_count, 
              " | Neg: ", neg_count
            )
          } else {
            direction_matrix[i, j] <- "Pos: 0 | Neg: 0"
          }
        } else {
          # Count taxa with same direction between two different treatments
          pos_i <- direction_data$Positive_Taxa[[i]]
          neg_i <- direction_data$Negative_Taxa[[i]]
          pos_j <- direction_data$Positive_Taxa[[j]]
          neg_j <- direction_data$Negative_Taxa[[j]]
          
          # Find taxa that are positive in both treatments
          same_pos <- length(intersect(pos_i, pos_j))
          # Find taxa that are negative in both treatments  
          same_neg <- length(intersect(neg_i, neg_j))
          
          direction_matrix[i, j] <- paste0(
            "Same Pos: ", same_pos, 
            " | Same Neg: ", same_neg
          )
        }
      }
    }
    
    direction_df <- as.data.frame(direction_matrix) %>%
      tibble::rownames_to_column("Treatment") %>%
      dplyr::arrange(factor(Treatment, levels = treatment_order))
    
    direction_table <- direction_df %>%
      gt::gt() %>%
      gt::tab_header(
        title = "Direction of Effect Comparison",
        subtitle = paste(section_title, "- Shared taxa with same direction of effect")
      ) %>%
      gt::tab_style(
        style = cell_text(weight = "bold"),
        locations = cells_column_labels()
      )
    
    # Apply treatment colors to column headers and row headers
    for (i in seq_along(treatments)) {
      trt <- treatments[i]
      if (trt %in% names(treatment_color_scale)) {
        trt_color <- treatment_color_scale[trt]
        
        # Color column header
        direction_table <- direction_table %>%
          gt::tab_style(
            style = cell_fill(color = trt_color),
            locations = cells_column_labels(columns = trt)
          ) %>%
          gt::tab_style(
            style = cell_text(color = "white", weight = "bold"),
            locations = cells_column_labels(columns = trt)
          )
        
        # Color row header (first column) for matching treatment
        direction_table <- direction_table %>%
          gt::tab_style(
            style = cell_fill(color = trt_color),
            locations = cells_body(
              columns = "Treatment",
              rows = Treatment == trt
            )
          ) %>%
          gt::tab_style(
            style = cell_text(color = "white", weight = "bold"),
            locations = cells_body(
              columns = "Treatment",
              rows = Treatment == trt
            )
          )
      }
    }
  } else {
    direction_table <- NULL
  }
  
  return(list(
    summary_table = summary_table,
    top10_table = top10_table,
    overlap_table = overlap_table,
    direction_table = direction_table
  ))
}

Import Data

ps.tmp <- readRDS("/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds")


ps.cleaned <-
    ps.tmp %>%
        ## Update Metadata
        ps_rename(Time = Timepoint) %>%
        microViz::ps_mutate(
            Treatment = case_when(
                Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
                Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
                Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
                Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
                Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
                Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
                Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
                Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
                TRUE ~ "Unknown"
            ), .after = "Pathogen"
        ) %>%
        microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
        microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
        microViz::ps_filter(Treatment != "Unknown") %>%
        microViz::ps_mutate(
            History = case_when(
                Antibiotics + Temperature == 0 ~ 0,
                Antibiotics + Temperature == 1 ~ 1,
                Antibiotics + Temperature == 2 ~ 2,
            ), .after = "Treatment"
        ) %>%
        
        ## Additional metadata updates, factorizing metadata
        microViz::ps_mutate(
        # Create treatment code
            treatment_code = case_when(
              Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
              Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
            ),
            # Create treatment group factor
            treatment_group = case_when(
              Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
              TRUE ~ "Control"
            ),
            # Convert to factor with appropriate levels
            treatment_group = factor(treatment_group, 
                                   levels = c("Control", "Parasite", 
                                              "Antibiotics", "Antibiotics_Parasite",
                                              "Temperature", "Temperature_Parasite",
                                            "Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
                                   ),
            treatment_code = factor(treatment_code, levels = treatment_order),
            # Create time point factor
            time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
            # Create pathogen status factor
            pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
                                   levels = c("Unexposed", "Exposed")),
            # Create sex factor
            sex = factor(Sex, levels = c("M", "F"))
            )  %>%
    microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
      microViz::ps_mutate(Exp_Type = case_when(
          Treatment %in% c("A- T- P-", "A- T- P+")  ~ "No prior stressor(s)",
          Treatment %in% c("A+ T- P-", "A+ T- P+")  ~ "Antibiotics",
          Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
          Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
      )) %>%
      microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
  # Fix names for taxonomic ranks not identified
  microViz::tax_fix(suffix_rank = "current", anon_unique = T, unknown = NA) %>% 
  # Filter for any samples that contain more than 5000 reads
  microViz::ps_filter(sample_sums(.) > 5000) %>%
  # Any taxa not found in at least 3 samples are removed
  microViz::tax_filter(min_prevalence = 3, undetected = 0) %>%
  # Remove any unwanted reads
  microViz::tax_select(c("Mitochondria", "Chloroplast", "Eukaryota"), deselect = TRUE) %>%
  microViz::tax_select(c("Bacteria, Phylum"), deselect = TRUE) 


# ps.cleaned %>% microViz::samdat_tbl()

Variables

## Diversity Methods

diversity.method <- list()
diversity.method[["alpha"]] <- c("shannon", 
                                 "inverse_simpson", 
                                 "observed")
diversity.method[["beta"]] <- c("bray", "canberra")

## Stats/Plotting Variables

alpha.stats <- list() # save all stat results here
alpha.plots <- list() # save plots here

beta.stats <- list() # save all stat results here
beta.plots <- list() # save plots here

diffAbnd.stats <- list() # save all stat results here
diffAbnd.plots <- list() # save plots here

00) All

MaAsLin2

60 DPE

Heatmap

Tables
Summary
Summary of Significant Taxa by Treatment
All Treatments at 60 DPE
Treatment Total Significant Positive Effect Negative Effect
A- T- P+ 44 22 22
A+ T- P- 32 10 22
A+ T- P+ 52 13 39
A- T+ P- 34 8 26
A- T+ P+ 50 23 27
A+ T+ P- 31 21 10
A+ T+ P+ 64 39 25
Top 10
Top 10 Most Significant Taxa by Treatment
All Treatments at 60 DPE - Ranked by q-value
Taxon Treatment Coefficient Phylum q-value
Rhodovarius A- T- P+ 1.514 Pseudomonadota 3.235 × 10−9
X67 14_Genus A- T- P+ 3.838 NA 3.235 × 10−9
Polymorphobacter A- T- P+ 2.812 Pseudomonadota 8.272 × 10−8
Plesiomonas A- T- P+ −2.715 Pseudomonadota 1.745 × 10−7
Agromyces A- T- P+ −3.414 Actinomycetota 6.097 × 10−7
Culicoidibacter A- T- P+ −3.504 Bacillota 1.513 × 10−6
Vogesella A- T- P+ 3.508 Pseudomonadota 2.155 × 10−6
Pseudomonas A- T- P+ 1.887 Pseudomonadota 5.315 × 10−6
Rhodocyclaceae_Genus A- T- P+ 2.383 NA 8.385 × 10−6
Vampirovibrionaceae_Genus A- T- P+ 3.092 NA 9.089 × 10−6
Cloacibacterium A+ T- P- 3.222 Bacteroidota 2.930 × 10−9
Caulobacter A+ T- P- 3.099 Pseudomonadota 4.057 × 10−9
Gemmataceae_Genus A+ T- P- 4.521 NA 6.193 × 10−8
Barnesiellaceae_Genus A+ T- P- 2.655 NA 7.371 × 10−7
Agromyces A+ T- P- −3.381 Actinomycetota 2.833 × 10−6
Polyangiia_Genus A+ T- P- −1.087 NA 5.742 × 10−6
Eremiobacteria_Genus A+ T- P- −1.208 NA 4.638 × 10−4
Rhodospirillales_Genus A+ T- P- −0.365 NA 6.675 × 10−4
Bacteroidia_Genus A+ T- P- −3.282 NA 6.923 × 10−4
Piscinibacter A+ T- P- −0.696 Pseudomonadota 8.677 × 10−4
Acinetobacter A+ T- P+ 5.204 Pseudomonadota 1.602 × 10−11
Tundrisphaera A+ T- P+ −5.278 Planctomycetota 5.902 × 10−9
Halomonas A+ T- P+ 2.321 Pseudomonadota 6.097 × 10−7
Uliginosibacterium A+ T- P+ 1.617 Pseudomonadota 7.610 × 10−7
Polyangiia_Genus A+ T- P+ −1.087 NA 2.363 × 10−6
Paracoccaceae_Genus A+ T- P+ −2.928 NA 2.622 × 10−6
Roseomonas A+ T- P+ −3.359 Pseudomonadota 2.885 × 10−6
Vibrio A+ T- P+ 3.549 Pseudomonadota 8.277 × 10−6
Legionella A+ T- P+ −3.469 Pseudomonadota 9.368 × 10−6
Gemmobacter A+ T- P+ −2.481 Pseudomonadota 2.461 × 10−5
Culicoidibacter A- T+ P- −5.692 Bacillota 1.214 × 10−12
Agromyces A- T+ P- −3.753 Actinomycetota 6.097 × 10−7
Flavobacterium A- T+ P- −3.608 Bacteroidota 2.476 × 10−6
Legionella A- T+ P- −3.757 Pseudomonadota 9.423 × 10−6
Polyangiia_Genus A- T+ P- −1.087 NA 1.311 × 10−5
Ensifer A- T+ P- −3.281 Pseudomonadota 3.027 × 10−5
Aquihabitans A- T+ P- −3.249 Actinomycetota 7.778 × 10−5
Paenirhodobacter A- T+ P- −2.895 Pseudomonadota 1.157 × 10−4
Dinghuibacter A- T+ P- −3.151 Bacteroidota 1.323 × 10−4
Sutterellaceae_Genus A- T+ P- −3.665 NA 1.616 × 10−4
Vogesella A- T+ P+ 5.733 Pseudomonadota 1.182 × 10−13
Unknown_Family_Genus A- T+ P+ 1.662 NA 2.150 × 10−7
Plesiomonas A- T+ P+ −2.666 Pseudomonadota 7.610 × 10−7
Pseudomonas A- T+ P+ 2.089 Pseudomonadota 1.122 × 10−6
Culicoidibacter A- T+ P+ −3.524 Bacillota 3.117 × 10−6
Polyangiia_Genus A- T+ P+ −1.087 NA 3.900 × 10−6
Agromyces A- T+ P+ −3.260 Actinomycetota 4.394 × 10−6
Rubrivivax A- T+ P+ 3.068 Pseudomonadota 8.593 × 10−6
Polymorphobacter A- T+ P+ 2.358 Pseudomonadota 1.590 × 10−5
Paenirhodobacter A- T+ P+ −2.864 Pseudomonadota 4.843 × 10−5
Dongia A+ T+ P- 0.482 Pseudomonadota 4.733 × 10−9
Flavobacterium A+ T+ P- −5.032 Bacteroidota 3.755 × 10−8
Candidatus_Nomurabacteria_Genus A+ T+ P- 0.396 NA 1.133 × 10−5
Neochlamydia A+ T+ P- 2.532 Chlamydiota 1.244 × 10−4
Polyangiia_Genus A+ T+ P- −1.087 NA 3.883 × 10−4
Vampirovibrionaceae_Genus A+ T+ P- 3.326 NA 4.156 × 10−4
RBG 13 54 9_Genus A+ T+ P- 0.538 NA 7.651 × 10−4
Bosea A+ T+ P- −2.954 Pseudomonadota 9.486 × 10−4
Agromyces A+ T+ P- −3.011 Actinomycetota 1.271 × 10−3
Thermoactinomycetaceae_Genus A+ T+ P- 0.214 NA 2.018 × 10−3
Pararhodobacter A+ T+ P+ 0.597 Pseudomonadota 4.872 × 10−10
Polymorphobacter A+ T+ P+ 3.167 Pseudomonadota 6.193 × 10−8
Pseudoxanthomonas A+ T+ P+ 3.005 Pseudomonadota 8.272 × 10−8
Vampirovibrionaceae_Genus A+ T+ P+ 4.058 NA 1.639 × 10−7
Chthoniobacter A+ T+ P+ 2.782 Verrucomicrobiota 3.199 × 10−7
Flavihumibacter A+ T+ P+ 2.663 Bacteroidota 6.824 × 10−7
Candidatus_Megaira A+ T+ P+ 2.687 NA 1.122 × 10−6
Acetivibrio A+ T+ P+ 0.502 Bacillota 5.069 × 10−6
Rhodocyclaceae_Genus A+ T+ P+ 2.681 NA 6.610 × 10−6
Pseudomonas A+ T+ P+ 2.064 Pseudomonadota 7.831 × 10−6
Overlap
Taxa Overlap Between Treatments
All Treatments at 60 DPE - Number of shared significant taxa
Treatment A- T+ P+ A- T+ P- A+ T- P+ A+ T+ P+ A+ T- P- A- T- P+ A+ T+ P-
A- T- P+ 27 18 27 26 12 44 8
A+ T- P- 14 10 15 16 32 12 10
A+ T- P+ 28 21 52 27 15 27 9
A- T+ P- 14 34 21 17 10 18 8
A- T+ P+ 50 14 28 33 14 27 12
A+ T+ P- 12 8 9 12 10 8 31
A+ T+ P+ 33 17 27 64 16 26 12
Direction
Direction of Effect Comparison
All Treatments at 60 DPE - Shared taxa with same direction of effect
Treatment A- T+ P+ A- T+ P- A+ T- P+ A+ T+ P+ A+ T- P- A- T- P+ A+ T+ P-
A- T- P+ Same Pos: 1 | Same Neg: 14 Same Pos: 2 | Same Neg: 6 Same Pos: 2 | Same Neg: 18 Same Pos: 1 | Same Neg: 7 Same Pos: 1 | Same Neg: 12 Pos: 22 | Neg: 22 Same Pos: 3 | Same Neg: 13
A+ T- P- Same Pos: 15 | Same Neg: 18 Same Pos: 2 | Same Neg: 9 Same Pos: 7 | Same Neg: 21 Same Pos: 1 | Same Neg: 11 Pos: 10 | Neg: 22 Same Pos: 1 | Same Neg: 12 Same Pos: 9 | Same Neg: 18
A+ T- P+ Same Pos: 7 | Same Neg: 20 Same Pos: 0 | Same Neg: 8 Pos: 13 | Neg: 39 Same Pos: 0 | Same Neg: 12 Same Pos: 7 | Same Neg: 21 Same Pos: 2 | Same Neg: 18 Same Pos: 7 | Same Neg: 20
A- T+ P- Same Pos: 3 | Same Neg: 8 Pos: 8 | Neg: 26 Same Pos: 0 | Same Neg: 8 Same Pos: 2 | Same Neg: 7 Same Pos: 2 | Same Neg: 9 Same Pos: 2 | Same Neg: 6 Same Pos: 2 | Same Neg: 5
A- T+ P+ Pos: 23 | Neg: 27 Same Pos: 3 | Same Neg: 8 Same Pos: 7 | Same Neg: 20 Same Pos: 3 | Same Neg: 11 Same Pos: 15 | Same Neg: 18 Same Pos: 1 | Same Neg: 14 Same Pos: 10 | Same Neg: 16
A+ T+ P- Same Pos: 10 | Same Neg: 16 Same Pos: 2 | Same Neg: 5 Same Pos: 7 | Same Neg: 20 Same Pos: 1 | Same Neg: 8 Same Pos: 9 | Same Neg: 18 Same Pos: 3 | Same Neg: 13 Pos: 21 | Neg: 10
A+ T+ P+ Same Pos: 3 | Same Neg: 11 Same Pos: 2 | Same Neg: 7 Same Pos: 0 | Same Neg: 12 Pos: 39 | Neg: 25 Same Pos: 1 | Same Neg: 11 Same Pos: 1 | Same Neg: 7 Same Pos: 1 | Same Neg: 8

01) What is the host-microbiome response to parasite exposure?

MaAsLin2

60 DPE

Heatmap

Tables
Summary
Summary of Significant Taxa by Treatment
Parasite exposure only at 60 DPE
Treatment Total Significant Positive Effect Negative Effect
A- T- P+ 30 16 14
Top 10
Top 10 Most Significant Taxa by Treatment
Parasite exposure only at 60 DPE - Ranked by q-value
Taxon Treatment Coefficient Phylum q-value
Plesiomonas A- T- P+ −2.689 Pseudomonadota 1.365 × 10−9
Culicoidibacter A- T- P+ −3.484 Bacillota 1.339 × 10−7
X67 14_Genus A- T- P+ 3.838 NA 6.540 × 10−7
Vogesella A- T- P+ 3.508 Pseudomonadota 1.034 × 10−5
Polymorphobacter A- T- P+ 2.457 Pseudomonadota 3.985 × 10−5
Rhodocyclaceae_Genus A- T- P+ 2.383 NA 1.010 × 10−4
Pseudomonas A- T- P+ 1.887 Pseudomonadota 1.035 × 10−4
Paenirhodobacter A- T- P+ −2.755 Pseudomonadota 1.725 × 10−4
Candidatus_Megaira A- T- P+ 2.076 NA 1.725 × 10−4
Gemmobacter A- T- P+ −2.100 Pseudomonadota 3.325 × 10−4
Overlap
## NULL
Direction
## NULL

02) Is the host-microbiome response to parasite exposure historically contingent?

MaAsLin2

60 DPE

Heatmap

Tables
Summary
Summary of Significant Taxa by Treatment
Parasite Exposure with Prior Stressors at 60 DPE
Treatment Total Significant Positive Effect Negative Effect
A+ T- P+ 23 4 19
A- T+ P+ 3 1 2
A+ T+ P+ 15 14 1
Top 10
Top 10 Most Significant Taxa by Treatment
Parasite Exposure with Prior Stressors at 60 DPE - Ranked by q-value
Taxon Treatment Coefficient Phylum q-value
Gemmataceae_Genus A+ T- P+ −3.424 NA 1.197 × 10−4
Rhodovarius A+ T- P+ −1.392 Pseudomonadota 4.207 × 10−4
Mycobacterium A+ T- P+ −3.154 Actinomycetota 4.207 × 10−4
X67 14_Genus A+ T- P+ −2.658 NA 1.052 × 10−3
Uliginosibacterium A+ T- P+ 1.617 Pseudomonadota 1.512 × 10−3
Vampirovibrionaceae_Genus A+ T- P+ −2.845 NA 1.552 × 10−3
Hyphomicrobiales_Incertae_Sedis_Genus A+ T- P+ −3.329 NA 1.928 × 10−3
Rubrivivax A+ T- P+ −2.481 Pseudomonadota 4.265 × 10−3
Roseomonas A+ T- P+ −2.355 Pseudomonadota 4.655 × 10−3
Pseudobdellovibrionaceae_Genus A+ T- P+ −2.549 NA 5.035 × 10−3
X67 14_Genus A- T+ P+ −3.179 NA 1.147 × 10−4
Unknown_Family_Genus A- T+ P+ 1.662 NA 1.197 × 10−4
Rhodovarius A- T+ P+ −1.146 Pseudomonadota 6.600 × 10−3
Pararhodobacter A+ T+ P+ 0.597 Pseudomonadota 5.204 × 10−5
Obscuribacteraceae_Genus A+ T+ P+ 1.671 NA 5.204 × 10−5
Pseudoxanthomonas A+ T+ P+ 2.824 Pseudomonadota 1.147 × 10−4
Rhodococcus A+ T+ P+ 3.105 Actinomycetota 2.740 × 10−4
Rhodovarius A+ T+ P+ −1.514 Pseudomonadota 4.207 × 10−4
Chitinibacter A+ T+ P+ 1.787 Pseudomonadota 2.059 × 10−3
Acetivibrio A+ T+ P+ 0.335 Bacillota 4.644 × 10−3
Simkaniaceae_Genus A+ T+ P+ 0.785 NA 5.205 × 10−3
Candidatus_Obscuribacter A+ T+ P+ 2.227 NA 5.735 × 10−3
Flavihumibacter A+ T+ P+ 1.883 Bacteroidota 6.600 × 10−3
Overlap
Taxa Overlap Between Treatments
Parasite Exposure with Prior Stressors at 60 DPE - Number of shared significant taxa
Treatment A+ T+ P+ A- T+ P+ A+ T- P+
A+ T- P+ 2 2 23
A- T+ P+ 1 3 2
A+ T+ P+ 15 1 2
Direction
Direction of Effect Comparison
Parasite Exposure with Prior Stressors at 60 DPE - Shared taxa with same direction of effect
Treatment A+ T+ P+ A- T+ P+ A+ T- P+
A+ T- P+ Same Pos: 0 | Same Neg: 1 Same Pos: 0 | Same Neg: 2 Pos: 4 | Neg: 19
A- T+ P+ Same Pos: 1 | Same Neg: 1 Pos: 1 | Neg: 2 Same Pos: 0 | Same Neg: 2
A+ T+ P+ Pos: 14 | Neg: 1 Same Pos: 1 | Same Neg: 1 Same Pos: 0 | Same Neg: 1

03) Is host-microbiome recovery historically contingent?

MaAsLin2

60 DPE

Heatmap

Tables
Summary
Summary of Significant Taxa by Treatment
No Parasite Exposure with Prior Stressors at 60 DPE
Treatment Total Significant Positive Effect Negative Effect
A+ T- P- 27 9 18
A- T+ P- 23 3 20
A+ T+ P- 23 17 6
Top 10
Top 10 Most Significant Taxa by Treatment
No Parasite Exposure with Prior Stressors at 60 DPE - Ranked by q-value
Taxon Treatment Coefficient Phylum q-value
Gemmataceae_Genus A+ T- P- 4.521 NA 6.162 × 10−8
Caulobacter A+ T- P- 2.884 Pseudomonadota 2.900 × 10−5
Cloacibacterium A+ T- P- 2.897 Bacteroidota 1.634 × 10−4
Agromyces A+ T- P- −3.210 Actinomycetota 4.972 × 10−4
Flavobacterium A+ T- P- −2.537 Bacteroidota 5.166 × 10−4
Barnesiellaceae_Genus A+ T- P- 2.655 NA 9.999 × 10−4
Bacteroidia_Genus A+ T- P- −2.927 NA 1.141 × 10−3
Cytophagales_Genus A+ T- P- −3.121 NA 2.136 × 10−3
Polyangiia_Genus A+ T- P- −1.087 NA 4.553 × 10−3
Pseudorhodoplanes A+ T- P- 1.316 Pseudomonadota 9.177 × 10−3
Culicoidibacter A- T+ P- −5.692 Bacillota 1.072 × 10−11
Flavobacterium A- T+ P- −3.608 Bacteroidota 1.238 × 10−6
Legionella A- T+ P- −3.615 Pseudomonadota 2.823 × 10−5
Ensifer A- T+ P- −3.026 Pseudomonadota 4.583 × 10−5
Acinetobacter A- T+ P- 2.971 Pseudomonadota 5.118 × 10−5
Aquihabitans A- T+ P- −3.198 Actinomycetota 6.435 × 10−5
Agromyces A- T+ P- −3.558 Actinomycetota 2.082 × 10−4
Paenirhodobacter A- T+ P- −2.722 Pseudomonadota 2.698 × 10−4
Dinghuibacter A- T+ P- −2.446 Bacteroidota 7.574 × 10−4
Sutterellaceae_Genus A- T+ P- −3.057 NA 1.224 × 10−3
Flavobacterium A+ T+ P- −5.032 Bacteroidota 1.431 × 10−8
Vampirovibrionaceae_Genus A+ T+ P- 2.869 NA 2.348 × 10−4
Dongia A+ T+ P- 0.482 Pseudomonadota 4.143 × 10−4
Bosea A+ T+ P- −2.717 Pseudomonadota 4.632 × 10−4
alphaI_cluster A+ T+ P- 1.923 NA 6.001 × 10−3
Neochlamydia A+ T+ P- 2.532 Chlamydiota 6.481 × 10−3
Chlamydiales_Genus A+ T+ P- 2.443 NA 8.308 × 10−3
Crenobacter A+ T+ P- −2.192 Pseudomonadota 9.177 × 10−3
Candidatus_Nomurabacteria_Genus A+ T+ P- 0.372 NA 1.089 × 10−2
Gaiellales_Genus A+ T+ P- 0.676 NA 1.694 × 10−2
Overlap
Taxa Overlap Between Treatments
No Parasite Exposure with Prior Stressors at 60 DPE - Number of shared significant taxa
Treatment A- T+ P- A+ T+ P- A+ T- P-
A+ T- P- 7 5 27
A- T+ P- 23 3 7
A+ T+ P- 3 23 5
Direction
Direction of Effect Comparison
No Parasite Exposure with Prior Stressors at 60 DPE - Shared taxa with same direction of effect
Treatment A- T+ P- A+ T+ P- A+ T- P-
A+ T- P- Same Pos: 0 | Same Neg: 3 Same Pos: 0 | Same Neg: 5 Pos: 9 | Neg: 18
A- T+ P- Pos: 3 | Neg: 20 Same Pos: 1 | Same Neg: 4 Same Pos: 0 | Same Neg: 3
A+ T+ P- Same Pos: 1 | Same Neg: 4 Pos: 17 | Neg: 6 Same Pos: 0 | Same Neg: 5